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Abstract We present a novel numerical method that allows the calculation of 
nonlinear force-free magnctostatic solutions above a boundary surface on which 
only the distribution of the normal magnetic field component is given. The 
method relies on the theory of force-free electrodynamics and applies directly 
to the reconstruction of the solar coronal magnetic field for a given distribution 
of the photospheric radial field component. The method works as follows: we 
start with any initial magnctostatic global field configuration (e.g. zero, dipole), 
and along the boundary surface we create an evolving distribution of tangential 
(horizontal) electric fields that, via Faraday's equation, give rise to a respective 
normal field distribution approaching asymptotically the target distribution. At 
the same time, these electric fields are used as boundary condition to numerically 
evolve the resulting electromagnetic field above the boundary surface, modelled 
as a thin ideal plasma with non-reflecting, perfectly absorbing outer boundaries. 
The simulation relaxes to a nonlinear force-free configuration that satisfies the 
given normal field distribution on the boundary. This is different from existing 
methods relying on a fixed boundary condition - the boundary evolves toward 
the a priori given one, at the same time evolving the three-dimensional field 
solution above it. Moreover, this is the first time a nonlinear force- free solution is 
reached by using only the normal field component on the boundary. This solution 
is not unique, but depends on the initial magnetic field configuration and on the 
evolutionary course along the boundary surface. To our knowledge, this is the 
first time that the formalism of force-free electrodynamics, used very successfully 
in other astrophysical contexts, is applied to the global solar magnetic field. 
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1. Introduction 

Force- Free Electrodynamics (hereafter FFE) is a formal name for time-dependent 
electromagnetism in an ideal plasma with negligible inertia (p = 0) and negligible 
gas pressure (f3 = 0). The formalism of FFE has been developed for various 
rclativistic astrophysical applications (pulsars, astrophysical jets, gamma-ray 
bursts, etc.) where the plasma supports electric currents and electric charges 
QGruzinov, 1999[ ). The equations of FFE are Maxwell's equations with nonzero 
electric fields as follows: 

— = cV x B - AttJ , (1) 

¥ = - cVxE > (2) 

V • B = . (3) 
These equations are coupled by the ideal MHD condition 

E ■ B = , (4) 

implying that in an ideal plasma E and B are everywhere perpendicular, and 
the force-free condition 

p e E+ ij x B = . (5) 

c 

Here, J is the electric current density, and p e = (47r) _1 V-E is the electric charge 
density. IGruzinovl () 1999)) showed that it is possible to solve for J in the above 
set of equations and thus express the electric current density as a function of the 
electric and magnetic field, namely 

c„^ExB c(B-VxB-E-VxE)^ 
J = — VE +—- 5 -B. (6 

4tt B 2 4tt B 2 w 

One can then numerically integrate Maxwell's equations (eqs [TJ [2]) to obtain 
the time evolution of the electric and magnetic fields. We developed a three- 
dimensional (hereafter 3D) finite-difference time-domain (FDTD) cartesian FFE 
code with non-reflecting, perfectly absorbing outer boundaries applied success- 
fully to the 3D structure of the pulsar magnctosphcre (Kalapot harakos and ContopoulosH 
{MM see also lYeel fTMHl |Spitkovskyl . 

We then realized that the same formalism may be applied to follow the 
temporal evolution of any open force-free ideal-MHD system from an initial 
to a final magnetostatic configuration (i.e. with zero initial and final electric 
fields) when the magnetic field distribution at a boundary evolves for a finite 
number of time-iteration steps. After sufficient time, when the magnetic field at 
the boundary will have reached a given distribution, and all electric fields will 
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have relaxed to zero, the above expressions degenerate into the following limited 
set of steady-state FFE equations, 

J = ^-VxB, (7) 

V • B = , and (8) 

J x B = . (9) 

A solution of the system of eqs. (JT])-© satisfying a given normal-field distri- 
bution on the boundary can be viewed as the final FFE equilibrium solution. 
Given the ill-posed nature of this problem (only the normal field component on 
the boundary is known), this solution is not unique and depends on the initial 
magnetostatic configuration and on the course toward equilibriun^]. 

The above methodology can be applied to the inner solar corona, that is, 
between 1 and 2 — 2.5r Q (solar radii). In this environment, motions in the line- 
tying photosphere largely but nonlincarly dictate the evolution in the magnetic 
field lines above, along which plasma is assumed 'frozen in' and forced to follow 
their path (j3 1). At the edge of the inner corona, known as the 'source surface' 
(Schattc n, Wilcox, and Ness, 1969 ) the plasma again approaches equipartition 
with the magnetic field. Beyond the source surface the plasma starts dominat- 
ing again in the form of the solar wind. Therefore, between lr Q (photosphere) 
and ~ 2.5r (source surface) the magnetic field is thought to approximately 
satisfy eqs. (0-© along with the divergence-free condition, cq. ((3]), and the 
ideal MHD condition, eq. ((3]). In discussing the global solar field, and driven by 
the necessity to use the observed photospheric magnetic field as the boundary 
condition, we ignore the fact that the latter is, in fact, forced, rather than 
force- free QMetcalf et al, 1995| |Georgoulis and LaBonte, 2004[ ). 

Routine measurements of the global photospheric field provide only the line- 
of-sight magnetic field component i?z,os( r o, @i 4 1 ) f° r given heliocentric spherical 
coordinates (9, <^)q This is typically accomplished by synoptic (Carrington) 
maps of Blos where the heliographic latitude and the Carrington longitude 
can be transformed into 9 and (f>, respectively. Synoptic maps take one solar 
rotation (~ 27 days at the equator) to complete one Carrington rotation (0° 
- 360° in Carrington longitude) but are again assumed to capture a snap- 
shot of the photospheric magnetic field for practical reasons (see, however, the 
attempt to evolve the photospheric field by means of a flux dispersal model 
( |Schrijver and De Rosa, 2003[ ) and use the modified boundary for global field 
extrapolations) . 



1 The above statements are correct mathematically. In practice, however, as shown in § 3, 
numerical dissipation sets in and the system evolves adiabatically through a sequence of 
nonlinear force-free magnetostatic equilibria with decreasing magnetic free energy towards 
a unique, minimum-energy and hence current-free (potential) equilibrium. 

2 The Vector SpectroMagnetograph (VSM; |Henney et q?.|l2009ll of the Synoptic Optical Long 
Term Investigations of the Sun (SOLIS) facility and the Helioseismic and Magnetic Imager 
(HMI) onboard the Solar Dynamics Observatory (SDO) is and will be providing, respectively, 
the first full-disk photospheric vector magnetograms. 
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Determining the distribution of coronal magnetic fields and subsequent elec- 
tric currents is, per the above, a formidable task. Global solar field extrapolations 

are current-free (potential) in their majority (Altschuler and Newkirk, 1969 Wang a nd Sheeley, 1992| | 



Luhmann et al., 2002; Schrijver and De Rosa, 20031. Some magnetohydrodynamic| 



modelling efforts are also present in the literature ( Linker et al., 1999]|Riley, Linker, and Mikic, 2001H 
|Roussev et al., 2003[ ) but they require immense computing resources. To our 
knowledge, the only technique providing a global solar nonlinear force-free field is 
the optimization method of T. Wiegelmann and collaborators. A magnctostatic 
version of this method was proposed by Wiegelmann (2007). The main idea was 
to simultaneously minimize the Lorentz force (V x B) x B and the divergence (V- 
B) of the field in the extrapolation volume in spherical coordinates - a generaliza- 



tion of the Cartesian implementation of Wiegelmann (2004). Efforts to further 



refine the method are discussed in |Tadesse, Wiegelmann, and Inhester] (|2009p 
while a further generalization into a magnetohydrostatic nonlinear force-free 
solution is described by | Wiegelmann et qL|(|2"0"0"T|) and lRuan et ctZ.I (|2008|) . These 
optimization efforts, however, generally require the vector magnetic field on the 
photospheric boundary. We hereby propose an alternative approach reaching a 
nonlinear force-free magnetostatic solution by only using the normal (radial) 
field component on the boundary, as follows: 

i) Start with a simple static potential magnetic field distribution such as B = 

everywhere, or a dipolar field of any direction, etc. 
ii) Evolve the vertical photospheric magnetic field by imposing a distribution 

of horizontal photospheric electric fields using Faraday's equation (eq. 0). 

This distribution is chosen such that the resulting photospheric magnetic 

field asymptotically (t oo) approaches a given photospheric magnctogram 

(see next section). 

Hi) During the evolution of the photospheric magnetic field, force-free electrody- 
namic waves are injected from the photosphere into the corona. Assuming 
non-reflecting, perfectly absorbing conditions at large distances, these waves 
will be absorbed by the outer boundaries and will not re-enter the computa- 
tional domain. The moving electric charges carried by these waves begin to 
establish a nonlinear network of coronal electric currents. 

iv) Gradually, as the photospheric magnetic field distribution approaches the 
target, photospheric electric fields will correspondingly decrease, and as a 
consequence magnetospheric clcctrodynamic waves, electric fields, and electric 
charges will also asymptotically fade. The nonlinear network of coronal electric 
currents, however, will survive. 

v) When electric fields vanish everywhere, we relax to a force- free configuration 
that closely matches the given photospheric boundary condition. 

At this initial stage we cannot claim that our solution is representative of the 
actual global coronal field. In fact, and as argued above, our solution depends on 
the initial condition and on the course followed to reproduce the given boundary 
condition (see, however, footnote 1). To achieve a physically meaningful equi- 
librium solution we need to apply additional physical arguments that will be 
the subject of a future publication. This work simply lays out the details of 
the methodology and presents a few characteristic test cases. Notice 1 that our 
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methodology does not apply in studying the dynamical evolution of the solar 
corona between successive photospheric configurations (the characteristic speed 
of propagation of information in the solar corona, the Alfven speed, is on the 
order of 1,000 km/sec, whereas the characteristic speed of propagation of infor- 
mation in our electrodynamic 'corona', the speed of light, is ~ 300 times faster). 
It does apply, though, in determining the equilibrium coronal configuration. 

Our method has certain similarities with 'stress-and-relax' schemes proposed 
by e.g. lRoumeliot is (1996), Valori, Klie m, and Keppens (|2005p . in that we both 
act on the boundary magnetic field distribution, and then evolve the magne- 
tosphcric field towards stationarity. Our respective technical implementations, 
though, arc different (e.g. they work with the vector potential, we work only 
with the magnetic field; they 'stress' the horizontal boundary magnetic field 
component, we evolve the radial one; they rely on viscous damping in order to 
settle to stationarity, we implement absorbing outer boundary conditions that 
remove any non-stationarity from our computational box). 

In § 2 we present our numerical implementation of a photospheric electric field 
distribution which guarantees that the distribution of the radial photospheric 
magnetic field component will asymptotically approach that of a given solar 
magnetogram. In § 3 we follow the electrodynamic response of the corona under 
the action of the above photospheric electric field distribution, and present some 
representative solutions for various initial conditions. We also present the various 
quantities used to monitor the course toward equilibrium. Our conclusions are 
summarized and discussed in S 4. 



2. The photospheric magnetic field 

Let us first assume a local 2D staggered cartesian grid (x, y) of size element 
5 in the photosphere, with normal magnetic field B z defined at the center of 
each cell and tangential electric fields E x ,E y defined on the corresponding cells' 
edges, as shown in Fig. [T^,. Each cell is characterized by a vector position 
In such staggered mesh, total magnetic flux is identically conserved when we 
evolve magnetic fields through Faraday's equation. 

At t = 0, we populate our 2D grid with a particular force-free magnetostatic 
field configuration B z (i,j;t = 0) that can be anything (i.e., zero, dipole-like 
of any direction and strength, quadrupole-like, etc.). Our aim is to impose a 
particular horizontal electric field distribution that will asymptotically evolve 
B z toward a target magnetic field configuration Br(i,j)- We chose the following 
procedure: given the photospheric magnetic field distribution B z (i,j]t) at each 
time step t, we scan the full photospheric grid and at each cell position we 

add four electric field components at the four edges of the cell such that 

E x (i,j;t) -> E x (i,j;t)-f[B T (i,j)-B z (iJ;t)], 

E y (i,j;t) -> E y (i,j;t) + f[B T (i,j)-B z (i,j;t)}, 

E x (i,j + l;t) -> E x (i,j + l;t) + f[B T (i,j)-B z (i,j;t)}, 

E y (i+l,j;t) -> E y (i + l,j;t)-f[B T (i,j)-B z (i,j;t)], (10) 
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Figure 1. Implementation of electric field update algorithm (eqs. HOC at t = in a 2D 
staggered cartesian grid. Electric fields arc defined along grid edges. Magnetic fields arc defined 
in cell centers. Nonzero equal-magnitude target magnetic field Bi>(i,j) only at positions shown 
with circles. Circles with dots: Bj- outwards from the plane. Circles with crosses: Bj- into the 
plane. The grid is scanned left to right, top to bottom in frames a-f respectively. 



as shown in Fig. [TJl. Here / is a free positive numerical factor used to adjust 
the rate of convergence. We chose / = 0.5, as convergence becomes numerically 
unstable for / > 0.5 and too slow for / < 0.5. At each time step t, we start with 
E x (i,j;t) = E y (i,j;t) = everywhere, and then scan the 2D grid and update 
electric fields according to eqs. (|TU1) . If we only had one computational grid-cell, 
the z component of Faraday's equation (eq. [2J would be written in a first-order 
discretized form as follows: 

dB z (t) Acf[B T -B z (t)] 

dt 6 1 ' 

This equation would be integrated numerically to obtain the asymptotic solution 
B z (t —> oo) = Bt- When the algorithm is applied sequentially to the full hori- 
zontal grid, though, electric fields on edges that correspond to neighboring cells 
will be updated twice as the algorithm is applied in both cells. In the sketch of 
Fig. [lb-f, one can see how the electric fields of neighboring cells are updated for 
a simple magnetic field distribution consisting of three plus three neighboring 
cells with equal and opposite magnetic fields. An important point to make here 
is that the process can be generalized in spherical coordinates (a step to take in 
the future) besides the Cartesian implementation shown here. 
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The above 2D procedure can be directly generalized in 3D. We use a cubic 
cartesian grid in order to avoid numerical artifacts associated with the axis 
of a spherical or cylindrical coordinate system. Note that in such a grid, the 
photosphere nowhere coincides with any of the grid cell phases. Wc define as our 
solar boundary surface the outer surface of the largest cartesian volume lying 
entirely within r Q . Given a target photospheric radial magnetic field distribution 
BT(rQ,9,4>)r in heliocentric spherical coordinates (see the discussion in § 4 on 
how Bt is obtained from a Carrington magnctogram) , we define the target 
magnetic field at the center of every boundary cell surface at positions (x, y, z) as 
the projection perpendicular to that surface of the field Bt{tq, 9, 0)(fo/r) 2 f = 
(Btx, Bt v , Btz) (for practical purposes, it is assumed here that the magnetic 
field inside a certain depth below the photosphere is purely radial, since our 
cartesian grid does not trace too accurately the solar surface). As before, we 
then scan the boundary surface and update electric fields as 

i,3,k;t) - f(B Tz (i,j,k) - B z (i,j,k;t)) , 
i,j,k;t) + f(B Tz (i,j,k) - B z (i,j,k;t)) , 
i,j + l,k;t) + f(B Tz (i,j,k)-B z (i,j,k;t)) , 
i + l,j,k;t)- f{B Tz (i, .7, k) - B z (i, j, k;t)), 

k;t) ~ f(B Tx (i,j,k) ~ B x (i,j, k;t)) , 
i,3,k;t) + f(B Tx (i,j,k) - B x (i,j,k;t)) , 
i,j,k + l;t) + f(B Tx {i, j, k) - B x (i, j, k; t)) , 
i,j + l,k;t) - f(B Tx (i,j,k) - B x (i, j, k;t)) , 

hj,k;t) - f(B Ty (i,j,k) - B y (i,j,k;t)) , 
i,3,k;t) + f{B Ty (i,j,k) - B y (i,j,k;t)) , 

k + l;t) + f(B Ty (i,j, k) - B y (i,j, k; t)) , 
i + l,j, k: t) - f(B Ty {i,j, k) - B y (i,j, k; t)) , 

(12) 

for the corresponding boundary cells that face in the z, x, and y direction 
respectively As in 2D, each boundary electric field component is updated twice 
from its two neighboring boundary surface grid cells. 

To test our algorithm we have used a synoptic magnctogram mainly referring 
to Carrington rotation 2009 recorded by the Michelson-Doppler Imager (MDI; 
IScherrer et a/.[[TM5)) onboard the Solar and Heliospheric Observatory (SoHO). 
The choice was deliberate: Carrington rotation 2009 corresponds to the period 
of October-November 2003 (the so-called 'Halloween 2003' period) when the Sun 
exhibited unusually high eruptive activity and the solar active-region belt was 
populated by several very complex active regions (perhaps an active-region nest, 
as proposed by IZhou et all I2007P . Using a solar dipole aligned with the dipole 
component of the given magnetogram as the initial condition, the photospheric 
magnetic field evolution toward the target is depicted in Fig. [21 Here, the entire 



E x (i,j,k;t) -» E x 

E y (i,j,k;t) — » E y 

E x (i,j + l,k\t) — » E x 

E y (i + l,j,k;t) -> E y 

E y (i,j,k;t) -> Ey 

E z (i,j,k;t) Eg 

E y (i,j, k+l;t) -> E y 

E z {i,j + l,k;t) -> E z 

E x (i,j,k;t) -> E x 

E z (i,j,k;t) -> E z 

E x (i,j,k + l;t) -> E x 

EJi + l,j,k;t) -> E z 
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solar disk is represented in (9, </>) heliocentric coordinates. Photospheric snap- 
shots are labelled by the number n of numerical time-integration steps, with 
the target magnetogram denoted with an infinite iteration number (n = oo). 
One may notice that the field evolution mimics flux emergence/submergence 
through the solar photosphere. Our procedure quickly reproduces the complex 
active regions present in the disk, but takes much longer to reproduce the weaker, 
large-scale field distribution in the polar regions. 

At this point, we should discuss the spatial resolution of our computational 
grid. Ideally we would like our uniform cubic grid to be able to capture both 
the details of the 'small-scale', strong active- region magnetic field and the large- 
scale structure of the much weaker coronal field. The first requirement calls 
for a size element 6 comparable to the resolution of the SoHO/MDI synoptic 
magnetogram, that is, <5 <~ 0.01x Q or 0.5° per pixel. Coupled with the second 
requirement that dictates a spatial extent of ~ 1.5rQ from the photosphere, it 
would give rise to a computational domain of ~ 10 8 cubic elements. For our 
calculation, that is performed on a typical desktop workstation, this grid size is 
infeasible. We thus reduced the original resolution of the synoptic magnetogram 
by a factor of 4 (8 ~ 0.04rQ or ~ 2° per pixel). If we were extrapolating for 
only a part of the solar disk, say, an active region, this restriction would be 
much less stringent, of course. As a result, the target synoptic magnetogram in 
Fig. [2]is shown with resolution degraded by a factor of 4. Although we cannot be 
certain about the degradation factors of other methods, say, the magnetostatic 
or magnctohydrostatic methods of T. Wicgelmann and colleagues, significant 
spatial degradation is, indeed, a necessity that one has to deal with when cal- 
culating the global solar magnetic field. The numerical time-integration step of 
our simulation is 0.5<5/c = 0.05 sec. Our simulations need about 24 hours on 
a desktop Intel i7 quad-core workstation to complete 10,000 time-integration 
steps. 

3. The coronal magnetic field 

Given the above boundary condition for the photospheric electric field we now 
follow the respective magnetic field evolution in the corona. This is treated as 
a pure electrodynamic problem that can be described through the formalism of 
FFE. We emphasize, though, that this is not a vacuum electrodynamic problem 
since electric currents and electric charges are allowed to develop everywhere in 
the corona to comply with the ideal MHD condition, eq. (@| . 

Electrodynamic waves will traverse the computational region informing the 
corona that the photosphere is changing. As they do so, they leave behind a 
modified coronal electric and magnetic field distribution. Similar waves are con- 
tinuously generated during photospheric magnetic flux emergence/submergence 
and travel out to larger heights in the corona to be absorbed in the non force-free 
region of the solar wind. This absorption is implemented numerically through a 
Perfectly Matched Layer (hereafter PML; |Berenger[ IM)4l ITM)l see also 
|Kalapotharakos and Contopoulos[ [2009[1 . that is, a non-reflecting, perfectly ab- 
sorbing boundary that mimics open space. 
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Figure 2. Evolution of the simulated photospheric normal field component in heliocentric 
(8, <j>) coordinates. Sub-plots labelled according to time-integration step. The same photo- 
spheric instances are shown with two different color scales. Left: linear color scale from -lkG 
to +lkG. Right: brighter areas: B r > 0; darker areas: B r < 0. Initial configuration: dipolar, 
with B r (9 = 0°) = —0.85 G. Target magnetogram labelled "oo". This time-sequence results 
in the equilibrium solution shown in Fig. \3\ 

As the distribution of the radial photospheric field approaches that implied 
by the target magnetogram, the photospheric electric fields will gradually fade. 
As a consequence, the coronal electrodynamic waves and electric fields will 
also fade, and the coronal magnetic field will gradually approach a force-free 
magnetostatic configuration satisfying the ideal MHD condition. In doing so, a 
nonlinear coronal current network develops. In Figs.JSJIH we show two examples 
of global coronal solutions reached for two different initial coronal field config- 
urations. In the lower right panels we plot the evolution with time-integration 
step of i) the total magnetic energy J B 2 dV over our computational volume V 
normalized to its asymptotic value obtained after n = 10, 000 time-integration 
steps (dashed line), ii) the average sine of the angle between the electric current 
J oc V x B and the magnetic field B weighted by the magnitudes of J and 
B (dotted line), and iii) the average absolute difference between the target 
and actual magnctograms normalized to the average magnitude of the target 
photospheric magnetic field (solid line). As can be seen, the first two monitors 
approach numerical convergence to within 10% in about 500 time-integration 
steps, but the third monitor takes roughly ten times as many steps to drop below 
the level of 10%. Indeed, the photospheric magnetic field needs longer integration 
times to settle to its asymptotic value, implying that surface electric fields take 



SOLA: Lowresolution.tex; 29 November 2010; 2:58; p. 9 



Contopoulos, Kalapotharakos, Georgoulis 




Figure 3. Global nonlinear force-free coronal configuration corresponding to the target 
photosphcric radial field of Carrington rotation 2009. Initial dipole polar magnetic field: 
B r (6 = 0°) = —0.85 G. Top panel: 3D configuration of magnetic field lines as viewed from the 
direction <j> = 0°. Field-line color: twist parameter a = (V X B)/B. Photosphcric color scale as 
in Fig. [2] Lower left panel: final configuration of magnetic fields projected onto a 2D meridional 
cut along <j> = 90°. Color scale: absolute value of the electric current J. Lower right panel: 
evolution of convergence monitors with numerical time-integration step. Dashed line: magnetic 
field energy normalized to its asymptotic value attained after 10,000 steps. Dotted line: average 
weighted sine of the angle between J and B. Solid line: average absolute difference between 
target and actual magnctograms over average absolute target magnetic field. Arrow: numerical 
integration step that corresponds to the plotted solution (this is when all convergence monitors 
fall below 10% of their asymptotic value). 
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Figure 5. Coronal configuration with nearly zero electric currents (potential) that corresponds 
to the same target photosphere as in figs. [3] & [J] obtained after 10,000 time integration steps. 



longer times to fade out. Currently we consider as equilibrium solution the state 
in which all convergence monitors fall within 10% of their asymptotic value. This 
practice will be reconsidered in future efforts. In the lower left panels we show the 
filial configuration of magnetic fields projected onto a 2D meridional cut along 
<fi = 90°, colored by the magnitude of the field-aligned electric current. In the 
top panels we show a 3D configuration of magnetic field lines as observed from 
the direction = 0°, colored according to the twist parameter a = (V x B)/B 
(magnetic field line color is redish where J is parallel to B, blucish where J is 
anti-parallel to B, and light green where J ~ 0). 

In Fig. [31 we begin with a dipole magnetic field counter-aligned with the solar 
axis with B r (9 = 0°) = —0.85 G, and aim for an actual global photospheric 
magnetic field obtained through part of Carrington magnctogram 2009. The 
chosen initial polar magnetic field value corresponds to the average one in the 
given magnetogram. By choosing this particular value we 'assist' the longer- 
term numerical convergence of our code. As can be seen in Fig. UJ the target 
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magnetogram is reproduced to within 10% in about 2,000 time integration steps 
(see also time sequence of Fig. [2]). This is when we plot the global coronal 
solution shown in Fig. [3J The mean weighted sine of the angle between the 
electric current and the magnetic field vectors is < 0.1, and the values of 
|V • B|/[(^) 2 + (^f ) 2 + (^r) 2 ]5 are on the order of 10~ 5 . The a parameter 
values are in the range ~ ±10 -3 Mm -1 (3 standard deviations away from its zero 
mean value) . Notice that such twist is not unreasonable in the solar corona. Here 
we should clarify that, although the photosphere has relaxed to its equilibrium 
configuration, continuing numerical integration beyond this stage results in a 
gradual disappearance of coronal electric currents due to the numerical diffusion 
inherent in any numerical code such as this one. In other words, coronal currents 
are continuously 'eroded' by numerical diffusivity. We know that, by the time 
we manage to reproduce our target photosphcric boundary condition to within 
10%, some amount of the coronal electric current will be lost due to numerical 
diffusivity. Beyond that point, the corona evolves through a sequence of mag- 
netostatic equilibria, each with smaller amount of electric current. After about 
10,000 time integration steps, we reach a configuration (Fig. [5]) with very small 
coronal electric currents (a ~ ±10 -4 Mm -1 ). This current-free configuration is 
uniquely determined for a closed volume or for a lower boundary condition and 
the semi-infinite space above it (e.g., |Aly] I1987[) . which is our approximation. 

Figure H] is similar to Fig. [31 only the initial magnetic field is dipolar counter- 
aligned with the solar axis with B r (6 = 0°) = —2 G. As before, we reproduce very 
quickly the main photosphcric active regions, but now the polar photospheric 
regions take longer to form, within about 4,000 time integration steps. One 
might expect that, due to the slower convergence, the resulting coronal electric 
currents are now slightly weaker than before (a values of the order ±0.75 x 
10 _3 Mm -1 ), and the two coronal configurations are somewhat different. If we 
continue the time integration for much longer times, coronal electric currents 
will also fade gradually, and the coronal configuration will approach the same 
potential configuration shown in Fig. [5] 

We also tested our algorithm in a much simpler configuration, namely that 
of a solar dipole without magnetic flux in the active-region belt. We used var- 
ious initial configurations (zero or dipolar). In all cases, the equilibrium solu- 
tion is almost indistinguishable from the vacuum magnctostatic dipole (poten- 
tial field), with very small field-aligned electric currents (a values of the order 
±10 _4 Mm -1 ). It seems that the absence of magnetic flux in the active-region 
belt gives rise to weak coronal electric currents and by the time the photosphere 
relaxes to its given boundary configuration numerical diffusivity has enough time 
to dissipate them effectively. In all cases, the average weighted value of the sine 
of the angle between the electric current and the magnetic field is < 0.1, and the 
maximum value of |V • B|/[(^) 2 + (^gf) 2 + (^) 2 }i is of the order 10 -5 . We 
obtained similar results when we considered different dipole orientations with 
respect to the axes of our numerical grid. This is evidence that our results are 
grid independent. 
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4. Discussion and conclusions 



Both solutions shown in Figs. |3] and 0] represent force-free magnetostatic con- 
figurations satisfying the ideal-MHD condition that correspond to the same 
magnetogram on the solar photosphere, therefore, they both represent solutions 
to our problem. The potential solution shown in Fig. [5] also represents a solution 
to our problem. Nature, however, chooses only one particular configuration. Our 
numerical approach suggests that our reconstructed field retains some memory 
of the prehistory of its evolution: it depends on the initial global configuration 
and on the course toward equilibrium. This is analogous to what we suspect 
happens in the real Sun as well, when different solar cycles succeed one an- 
other. If the target boundary does not have any well-formed active regions, 
i.e. like in a solar-minimum configuration, our simulation results in the unique 
potential-field solution between the photosphere and the source surface or for 
the semi-infinite space above the photosphere. This is also what full-fledged 
magnetohydrodynamic global models give, in this case approaching the results of 
potential-field models flRiley et al, 2006] ) . Obviously there is no straightforward 
answer to the question of how to choose one particular coronal field configuration 
from another. One way could be to introduce additional control monitors such as 
magnetic free energy in our computational volume or the Poynting flux through 
the photospheric boundary. Another way would be to make use of full vector 
magnetograms whenever those are available (for example, in eqs. (|12[) we could 
involve the three magnetic field components on the boundary toward the target 
components (Bt x , Bt v , Bt z ) of the observed field vector B; not only its radial 
component B r r). This major issue, that could fully constrain our solution, will 
be addressed in a forthcoming publication. 

Finally, we address a problem associated with the observational method via 
which the radial component of the photospheric magnetic field £?t(?"©, </0 is 
inferred from measurements of the line-of-sight magnetic field Blos{ t q>i 9, </0- 
We will only consider here the method applied to a Carrington synoptic mag- 
netogram, although our discussion below also applies to local observations of 
particular photospheric active regions, as well. In a Carrington magnetogram we 
arc observing the photosphere along solar meridional sections that contain our 
line of sight. As the Sun rotates, the full photosphere is covered and a global 
magnetogram is obtained under the assumption that the magnetic field structure 
remains unchanged during a full solar revolution (obviously, this is not quite 
true, but is the best one can do with present day satellite technology, avoiding 
also further interference with observational data). In order to obtain the radial 
magnetic field, a further assumption is made, namely that the magnetic field at 
the base of the corona is everywhere perpendicular to the photosphere, i.e., 

B T {r Q , &,(/)) = — . (13) 

sin 

This second assumption is also not totally accurate, but is the best one can do 
without knowledge of the photospheric field vector. The orientation of this field 
can be inferred from a full 3D reconstruction of the global coronal magnetic field, 
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and this is precisely what we are doing in the present work. We, therefore, pro- 
pose to include in a future publication one more step in the iteration procedure 
implemented along the photosphere, namely, 

i) At time t = 0, start with a purely radial target for the normal photospheric 
field as given by eq. ([T3f . 

ii) Apply the numerical integration steps and calculate the interim line-of-sight 
component i?LOS(ntim)( r 0j 4>\ t)- I n general, this will be different from 
B r {rQ,9,4>]t)sin9 because a nonzero horizontal magnetic field now develops 
in the photosphere. 

in) Rescale the target field Bt(t®, 6, <j>) by a factor reflecting the development of 
horizontal fields on the boundary, i.e., 



B T (r e ,0,< 



B L os{rQ,0,cj)) ( B r (r e ,6,(f>;t)smt 



li 'I III ) 

(r Q ,6>, </>;£) 



\-t>LOS(num){re,V, 9\t) , 



Here, the geometric term sin9, reflecting the assumption of purely radial 
field on the boundary in eq. (|13[) , is replaced by the more general term 

BLOS(num)/ B r . 

The above iterative procedure, implemented together with eqs. (|12p . contin- 
uously redefines the target radial photospheric field, Bt, gradually relaxing the 
assumption of purely radial fields in the photosphere. To our knowledge, this is 
the first time that the full 3D coronal field geometry has been taken into account 
in de-projecting the photospheric magnetic field. 

The quest toward choosing the most realistic global reconstruction of the solar 
coronal magnetic field at the lowest possible computational expense is, obviously, 
still ongoing. Despite the non-uniqueness of our solution we believe that the 
proposed method is interesting and may even turn out to become promising, 
due to the following two reasons: firstly, ours is the first technique that reaches a 
nonlinear force-free solution using only the radial photospheric field that is read- 
ily available from various sources. Attempting to match our force-free tangential 
evolution in the photosphere with the tangential component provided by SO- 
LIS/VSM or will be provided by SDO/HMI full-disk vector magnetograms may 
further constrain, or determine, hopefully, our solution. Secondly, our method 
alleviates the need for preprocessing the forced photospheric fields to better 
comply with a coronal force-free solution. Preprocessing introduces additional 



assumptions (e.g., Tadesse, Wicgelmann, and Inhester I2009j) so the fact that 



we asymptotically reach the target synoptic magnetogram at the same time 
ensuring a force-free solution above each interim boundary configuration may 
prove advantageous. Concerning our non-uniqueness problem, we will consider, 
classify, examine, and implement various revisions that may potentially present 
an opportunity to improve our global solar magnetic field reconstruction method. 

Acknowledgements This work uses synoptic solar magnetograms from SoHO/MDI. SoHO 
is a project of international cooperation between ESA and NASA. 
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